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The water molecule is a key ingredient in the formation of planetary systems, 
with the water snowline being a favourable location for the growth of 
massive planetary cores. Here we present Atacama Large Millimeter/ 
submillimeter Array data of the ringed protoplanetary disk orbiting the 
young star HL Tauri that show centrally peaked, bright emission arising 
from three distinct transitions of the main water isotopologue (H5°O). The 
spatially and spectrally resolved water content probes gas in a thermal range 
down to the water sublimation temperature. Our analysis implies a stringent 
lower limit of 3.7 Earth oceans of water vapour available within the inner 17 
astronomical units of the system. We show that our observations are limited 
to probing the water content in the atmosphere of the disk, due to the high 
dust column density and absorption, and indicate that the main water 
isotopologue is the best tracer to spatially resolve water vapour in 
protoplanetary disks. 


The water molecule is undoubtedly one of the most important mole- 
cular species in the whole Universe. Being an extremely efficient 
solvent, water had a key role in the emergence of life as we know it on 
our planet. For this reason, the chemical characterization of exo- 
planetary atmospheres is often focused on detecting this particular 
molecular species’ *. Formed by the common H and O atoms, water 
plays a fundamental role in the physics of the formation of planetary 
systems, due to its very high abundance in both gaseous and icy forms*”. 
Theoretical models predict that at the location of the phase transition 
from gaseous to solid form, dust grains can accumulate and grow very 
efficiently, promoting the fast formation of planetary cores. Across 
this particular radial location, called the ‘snowline’, grains can drasti- 
cally change their drift and fragmentation velocity, composition and 
opacity. In synergy with vapour radial diffusion‘, these physical discon- 
tinuities can lead to the accumulation and growth of dust grains into 


planetesimals”*. The position of the snowline also defines the chemistry 
of the available planet building blocks. Since the H,O molecule is the 
major elemental oxygen carrier in the disk, its desorption and freez- 
ing affect the elemental C/O ratio in both the gas and solid phases” ". 

Because of its large binding energy, the H,O transition from ice to 
gas happens a few astronomical units (au) from the young star where 
the midplane temperatures are in the range from 100 to 200 K, making 
it the last major ice component to sublimate. However, the proximity 
tothe host star makes the detection of the snowline complicated even 
in the closest star-forming regions. Both cold and warm water lines 
have been detected ina few disks by Herschel (see refs. 12,13 and refer- 
ences therein), Spitzer”, WST“ and ground-based observatories", 
but the low angular resolution did not allow robust inferences about 
the extent of the water snowline. Directly observing water emission 
from the ground is complicated by the high water vapour content of 
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Table 1| Observed H,O isotopologue transitions 


Transition Rest E, Channel r.m.s.° Flux Mask 
Frequency vo width radius” 
(GHz) (K) (kms") (mJybm") (mJykms”‘) (") 

p-H,O 183.31009 204.7 0.8 6.90 973+89 0.70 

313-220 

p-H20O 325.15290 469.9 1.0 3.71 1332+89° 0.70 

545-422 

o-H,O 321.22569 1,861.3 5.0 0.96 679+135 0.28 

1029-936 

p-H80 322.46517 4679 10 175 <33.64 0.70 

55-422 


The upper state energies are taken from ref. 51. The notation for the water energy levels in the 
vibrational ground state is Jç, x, “Obtained over one single channel. "Circular radius used to 
extract the line flux. ‘This flux measurement corrects for the absorption identified at ~1Okms". 
The flux obtained without accounting for absorption is 929+89 mJy kms”. “30 upper limit. 


the Earth’s atmosphere, resulting in strong telluric absorption. To 
circumvent this problem, most programmes at millimetre wave- 
lengths have focused on attempting the detection of the rarer H180 
and HDO isotopologues, leading to the clear detection of spatially 
resolved water isotopologue emission in the outbursting source V883 
Ori (ref. 19). In quiescent sources, except the candidate detections 
in the AS 205N and HL Tauri (HL Tau) disks, no detection of thermal 
emission at (submillimetre) wavelengths has been reported in 
the literature”. 

In this Article, we focus on the textbook case of HL Tau, the first 
protoplanetary disk imaged at very high angular resolution (~0.025”) 
with the Atacama Large Millimeter/submillimeter Array (ALMA)**. The 
disk shows a spectacular pattern of concentric rings. With the source 
being young (<1 Myr) and the dynamical stellar mass being relatively 
high (2.1 + 0.2 M,; ref. 25), the inner disk temperatures are expected to 
be warm, due to irradiation and accretion heating. Warm and hot water 
has been detected in HL Tau both by Herschel in the far infrared” and 
by ground-based high-resolution spectroscopy in the mid-infrared”, 
with the lines not being spatially resolved. With this Article, we report 
the detection of three rotational water lines in the inner regions of 
the HL Tau protoplanetary disk obtained with ALMA. The lines are 
spectrally resolved. Analysis of the interferometric data confirms that 
the extent of the water emission is confined within a prominent gap 
seen in the HL Tau continuum intensity. These data present spatially 
resolved images of the emission from the main water isotopologue 
(H,O) in a protoplanetary disk and pave the way to a new observa- 
tional strategy to characterize the water vapour content of terrestrial 
planet-forming regions. 


Observations 

We observed HL Tau in two different ALMA bands (band 5, originally 
developed with the goal of studying water in the local Universe”, and 
band 7) to target three transitions of water (two lines of para-water 
and one line of ortho-water): p-H,O 313-229 and 5,5-4», at 183.31 and 
325.15 GHz, respectively, and o-H,O 10,,-93, at 321.22 GHz. The first 
two lines are expected to trace warm water vapour outside the water 
snowline, while the third line is predicted to detect hot water within 
the water snowline*’*°. We also observed a rotational transition of 
the water isotopologue p-H0 at 322.46 GHz. The molecular coeffi- 
cients of the lines and sensitivity of the observations are reported in 
Tables land 2. 

After self-calibration, we imaged both the continuum and the 
continuum-subtracted lines with the CASA software”. The continuum 
images are shown in Fig. 1. All the H,O lines and the H80 line were 
imaged with CASA 6.2.1 with natural weighting, to maximize point 
source sensitivity. The 183 GHz line presents a synthesized beam of 
0.500” x 0.442” (position angle (PA) —3.0°) and was imaged witha 


channel width of 0.8 km ~. The resulting beam for the 325 GHz water 
line is 0.640” x 0.491” (PA -42.1°), with a 1 km s” channel spacing. For 
the high-excitation 321 GHz line, only one long baselines execution 
block was available, and the beam is 0.067” x 0.061” (PA 12.2°), witha 
5kms ‘channel spacing. The H1830 line was imaged with several differ- 
ent channel spacings. To have a one-to-one comparison with the main 
isotopologue line, in this Article we show the results with a channel 
width of 1 km s”. The resulting beam is 0.779” x 0.626” (PA -46.7°). 

The spectrum of the two lowest excitation lines was extracted 
over a circular area with 0.7” radius and is shown in Fig. 2. Both the 
183 and 325 GHz lines are clearly detected across multiple channels, 
with the lines being centred on the systemic velocity of the system 
(7.1kms7; refs. 32,33). The 325 GHz line shows an absorption compo- 
nent at~10 km s”, as seen for other lines at similar upper energy levels”. 
The 183 GHz line shows a width of -12 km s”, while the higher-frequency 
lines show a slightly broader emission. The 183 GHz line spectrum 
exhibits a peak signal-to-noise ratio (SNR) of ~9.6, with a root mean 
square (r.m.s.) of 13.2 mJy over the flux density extraction area. The 
325 GHz line spectrum instead shows a peak SNR of ~5.8, with ar.m.s. of 
14.0 mJy. The higher-energy line at 321 GHz does not showa detection 
when integrating over the same area, due tothe much smaller beam and 
the resulting high r.m.s. We thus extracted a spectrum over a circular 
area with radius 0.06”. The peak SNR is ~4.1, with ar.m.s. of 3.0 mJy. This 
transition shows the broadest line profile among the three detected 
lines, consistent with originating from the inner regions with higher 
Keplerian velocities. In all cases, the integrated intensity map shows 
a strong detection with a peak colocated with the dust continuum 
peak (Fig. 1). We obtain peak SNRs of 21.4, 19.8 and 8.1in the integrated 
intensity maps of the 183, 325 and 321 GHz lines, respectively. The line 
fluxes extracted over a circular area centred on the continuum peak 
are reported in Table 1. 

Theintensity-weighted velocity maps (moment 1 maps) of the 183 
and 325 GHz lines are shownin Fig. 1. Inthe high-SNR map of the 183 GHz 
line, disk rotation is clearly detected, with PA and systemic velocity in 
agreement with other brighter lines from the same dataset”, indicat- 
ing that a displacement of the photocenter of the blue-shifted and 
red-shifted channels is detected at the current resolution. 

For the H180 line in band 7, we extracted the spectrum with 
the same methodology described for the three main isotopologue 
lines, over a 0.7” radius circular area (Extended Data Fig. 1). A moment 
O map was computed over the same spectral range as the main isotopo- 
logue line. The line was not detected. 


The spatial distribution of water vapour 
The integrated intensity morphologies of the 183 and 321 GHz lines 
are remarkably different (Fig. 1), but so are the angular resolutions 
of the two moment maps. In order to derive radial profiles of the 
respective integrated intensities, we used two different approaches. 
First we focused on the highest-SNR line of the sample (at 183 GHz), 
which is also the coldest and therefore expected to trace the largest 
spatial extent. We averaged in frequency the interferometric data of 
the continuum-subtracted line in the same frequency range used to 
compute the moment O map. We then fitted the line integrated intensity 
visibilities with a simple Gaussian model using the galario package*™* 
after fixing the inclination and PA to 46.7 and 138.0°, respectively, 
as derived from high-angular-resolution continuum observations”. 
The fit converges well to a Gaussian with o, = 0.12 + 0.01” (Fig. 3 and 
Extended Data Fig. 2;0, is the standard deviation of the Gaussian func- 
tion). Ata distance of 140 pc, this corresponds too, = 16.8 + 1.4 au. For 
the 321 GHz line; the much higher angular resolution allowed us to 
compute the radial profile of the integrated intensity by azimuthally 
averaging the moment zero map after deprojecting it by the known 
inclination (using the GoFish package”). 

Figure 3 compares the two integrated intensity profiles with the 
band 7 azimuthally averaged continuum intensity profile. Both line 
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Table 2 | Observation IDs and execution block properties, including source integration time, median PWV column, 
bandpass, flux and phase calibrators and maximum baseline 


Program ID Lines Integration time (min) PWV (mm) Bandpass/flux calibrators Phase calibrators Maximum baseline (m) 
20171.01178.S p-H,0 343-220 46 0.2 J0423-0120 JO510+1800 1,398 

20171.01178.S p-H,0 545-422 31 0.4 J0538-4405 J0431+1731 3,637 

2017.1.01178.S 0-H,O 1045-936 33 0.5 J0519-4546 J0440+1437 8,547 

2022.1.00905.S p-H,O 545-422 100 0.3 J0423-0120 J0431+1731 500 
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Fig. 1| Continuum intensity and water vapour moment maps. a, Left, 1.7 mm shown due to low SNR. Right, zoom-in of continuum intensity, with [4,5,6,7,8]o 
continuum image of HL Tau. Centre, integrated intensity map ofthe 183 GHzwater contours of the 321 GHz line moment 0 map, with ø = 13.3 mJy beam’ km s”. The 
line. Right, intensity-weighted velocity map of the 183 GHz water line after 40 r.m.s. associated with the the integrated intensity maps of the 183 and 325 GHz 
clipping on individual channels, where disk rotation is clearly detected.b,Sameas lines are, respectively, 28.2 and 46.3 mJy beam“ km s~. Brightness intensity (K). 
a, for the 0.94 mm continuum and the 325 GHz water line. The intensity-weighted Integrated intensity (K km s”). Intensity weighted velocity along the line of sight 
velocity map in this case is computed after 30 clipping. c, left and centre, same (kms”). Distance from the phase center in Right Ascension (Delta RA). Distance 
as a, for the 0.94 mm continuum and the 321 GHz water line. No moment 1 map is fromthe same center in Declination (Delta Dec). 
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Fig. 2 | Water line spectra. a, Spectrum of the 183 GHz water line extracted over a 
circle with radius of 0.7” centred on the continuum peak. b, Spectrum of the 

325 GHz water line extracted over the same area, highlighting the mirrored 
(flipped) version of the spectrum with the dashed-dotted line. c, Spectrum of the 
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321 GHz water line extracted over a circle with radius of 0.06” centred on the 
continuum peak. The velocity range on the x axis is different in c. The grey dashed 
line in all panels shows the systemic velocity (Vsys) of 7.1 km s™. In the top left of 
each panel, 20 scale bars are reported for reference. 


profiles are clearly centrally peaked, indicating that the water vapour 
emission must originate above the optically thick continuum from the 
disk midplane. The lower excitation line is substantially more extended, 
showing that the water emission has a radially decreasing temperature 
(excitation) and that the high-excitation line is not optically thick out- 
side the central beam. The same line shows detectable emission out 
to 0.3” when boosting the SNR by azimuthally averaging. The 183 GHz 
temperature gradient is confirmed by fitting the high-SNR spectrum 
with a Keplerian disk model with a brightness temperature gradient, 
assuming that the line is close to being optically thick (Extended Data 
Fig. 3). Declining temperature profiles are preferred to flat ones, in 
agreement with the derived integrated intensity profile. 

The warm brightness temperature profile of the 183 GHz line, 
which far exceeds the midplane temperatures obtained by analysis of 
multiwavelength continuum data”, indicates that the water vapour 
we are tracing originates in the warm disk upper layers. This is further 
supported by the peak of the 321 GHz emission, whichis slightly shifted 
from the continuum peak (Fig. 1, bottom-right panel). Even though 
the two are consistent within the astrometric precision of the data, 
the apparent shift agrees with tracing water vapour on the side of the 
disk closer to the observer, which is in the northeast”, unocculted 
by the optically thick dust continuum. These findings set a stringent 
upper limit on the radius of the water snowline at 17 au. Amore accurate 


determination will require simultaneous forward-modelling of the 
radial and spectral profiles of all three (sub) millimetre lines with the 
aid of thermochemical codes. 


Water column density 

From the three main water lines, we computed the rotational diagram, 
under the assumption of optically thin emission and uniform excita- 
tion temperature across the energy range. While the very inner regions 
of the line emission are likely opaque, the bulk of the emission from 
the three lines cannot be optically thick, since this would imply a flux 
ratio that is equal to the square of the frequency ratio (in the Rayleigh- 
Jeans regime and local thermodynamical equilibrium), which we do 
not observe. While masing cannot be excluded as partly contributing 
to the emission, inthe high-spectral-resolution spectra of the 183 and 
325 GHzlines we have not identified individual narrow spectral features, 
which are in general a good signature for maser action together with 
high flux density. From the high densities of the inner disk, masing is 
not expected for the three lines analysed here, and the only contribu- 
tion could originate in the very upper layers at low volume densities 
where collisional quenching of the masing action is less probable. In 
the rotational diagram, we used degeneracy quantum numbers that 
account for a 3:1 ortho-to-para ratio, and a partition function that 
considers all states within the same assumption”. No rescaling of 
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Fig. 3 | Spatial extent of water emission. a, Integrated brightness temperature 
(T,) radial profile of the 321 Hz line, reconstructed integrated T, radial profile of 
the 183 GHz line and T, profile of the 0.94 mm continuum emission. The thick 
orange line shows the best-fit model of the 183 GHz line assuming a Gaussian 
integrated intensity profile. Thin lines show randomly sampled realizations of 
the posterior distribution. The lines in the top right show the beam major axis. 
For the 183 GHz line, it portrays the smallest spatial scales to which the uv plane 
analysis is sensitive, which are -2.3 smaller than the major axis of the respective 


uv distance (MA) 


natural beam. b, Real and imaginary part (Re and Im) of the recentred and 
deprojected visibilities (V) of the integrated intensity of the 183 GHz water line. 
The data points represent the weighted average of the visibilities within the 
60-kA-wide bins, with the scale bars showing the associated uncertainties 
(defined as UV È wi ), with w; being the weight of an individual visibility). The 
thick red line shows the best-fit model reproducing the profile shown ina. The 
units kA| and MA| express distance in units of the observed wavelength Al. 
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Fig. 4 | Rotational diagram. a, Rotational diagram of the three water lines, with line fluxes extracted as in the text. The fit does not lead to a unique solution, indicating 
that the assumption of uniform temperature is inadequate. b, Posterior distribution of the rotational diagram fit. The dashed lines indicate the 16th, 50th and 84th 


percentiles of the marginalized posterior distributions. 


the ortho- and para-line fluxes is needed to compute the total water 
column. We accounted for 10% absolute flux systematic uncertainties 
inthe line fluxes. 

The rotational diagram does not provide a unique solution for the 
rotation temperature, as shown in Fig. 4. The Monte Carlo Markov chain 
(MCMC) exploration individuates two distinct temperatures, which 
realistically indicate a continuous gradient in the excitation tempera- 
ture of the water vapour. Fitting either the two lower-energy lines 
or the two higher-energy lines separately, the rotational diagram indi- 
cates excitation temperatures of 214*5 and 789*}%” K, respectively, in 
line with the two classes of solutions obtained in the joint fit (Extended 
Data Fig. 4). The lower-temperature solution is sensitive to colder water 
vapour in the range of expected desorption temperatures of water ice 
in space, suggesting that the bulk of the water emission from the 
183 GHz line traces water gas in the proximity of the snow surface. The 
higher-temperature solution is driven by warm gas in the upper layers 
of the terrestrial planet-forming regions of the disk, which are well 
imaged by the high-resolution 321 GHz line integrated intensity map, 
and likely by the lines being close to be optically thick. The excitation 


temperature of the warm gas is in broad agreement with mid-infrared 
line fluxes in the 12.37-12.41 um range on the same source (in particular, 
the o-H,016,,;~15, ,, line with upper state energy Eup = 4,948 K (ref. 27)). 
These high temperatures allow for water vapour formed in situ via 
gas-phase reactions**””. 

Allthe rotational diagrams robustly constrain the column density 
of water gas (in the optically thin limit and above the optically thick 
continuum). The joint fit shows a total water column density 
logyo(Nehin/cm?) = 16.41*9 6° within 0.7” (-100 au) from the star. This 
corresponds to -3.7 Earth oceans (7.1 x 107 lunar masses) of 
water vapour. Given the optically thin assumption, this has to be 
considered a tight lower limit. Since most of the emission origi- 
nates from <0.12” (17 au), assuming that the entirety of it is confined 
within this radial range, we obtain an average column density 
of logy (Nenin/em?) ~ 18.10*0-0 , when accounting for the disk 
inclination. 

The non-detection of the H80 line can determine an upper 
limit of the average optical depth expected for the H,O 325 GHz line 
(7,0). Assuming that Th,o = 530 TyI80 (using the oxygen isotope ratio 
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Table 3 | Molecular coefficients of the observed H,O 
isotopologue transitions 


Transition gu Axy(s") Q Q Q Q 
(100K) (200K) (300K) (400K) 
p-H,O 3-2 7 359x10 351 97.4 1781 274.6 
PROS A 11 115x105 - p 2 ? 
o-H,O 10-9 63  6.21x10° - 2 = A 
1.05x10° - - - - 


p-H}8 O5s-42 1 


Values are from the Jet Propulsion Laboratory database, with radiative coefficients from 
ref. 52 and the updated partition function Q from the ExoMol database”, which well 
agrees with the one by ref. 52 in the temperature range explored in this Article. Only some 
representative values are reported in this table. The degeneracy quantum number of the 
ortho-state assumes an ortho-to-para ratio of 3. The partition function is computed over all 
states with the same ortho-to-para ratio. 


measured in the solar wind*°), we obtain that Tu o < 14. By turning 
the argument around, the optically thin assumption for the 325 GHz 
line provides an upper limit of N(H,O)/M(H}80) =40inHL Tau, well in 
agreement with oxygen fractionation levels in our own solar system”. 
The non-detection of H8Oshows that observational campaigns aimed 
at targeting water in inner disks of quiescent disks with ALMA should 
privilege the main isotopologue as a first choice, with follow-up 
observations of robust detections targeting HDO and H130 to derive 
more accurate column densities and set stringent limits on the 
water deuteration”. 

Assuming that the bulk of the water emission originates from 
<17 au from the star, we can compare the total mass of water (My, 0) to 
the mass of dust Maus estimated from multiwavelength continuum 
analysis of ALMA and Karl G. Jansky Very Large Array data in the same 
radial range*®. By using the dust surface density from this study, we 
obtain Maust = 13M garth, Which leads to a water-to-dust mass ratio of 
Muy,0/Mdust ~ 6 x 1075. This number is much lower than the expected 
water abundance in inner disks (water-to-dust mass ratio -10°). The 
optical thickness of water can only marginally alleviate the problem, 
given the non-detection of H380. Continuum subtraction could also 
marginally reduce the water brightness temperature” but in this case 
is not expected to reduce the water fluxes by more than a factor of 2 
(see Fig. 3). The low water-to-dust ratio further strengthens the 
interpretation that with ALMA, we are probing only the upper layers 
of the disk, above the optically thick screen of the dust continuum 
emission, which shows optical depths between ~5 and 10 at 0.9 mm 
within the inner 17 au (ref. 36). The large cavity seen in the HDO 
and H380 integrated intensity maps of V883 Ori (ref. 19) further 
supports that the observation of a large column of water is hindered 
by the optically thick continuum in the inner disk. 


Conclusions 

These new ALMA data reveal high-significance detections of three 
distinct rotational transitions of the main isotopologue of water 
vapour in the inner regions of the ringed HL Tau disk. These observa- 
tions pave the way to the characterization of the water content of the 
inner regions of protoplanetary disks. The tremendous angular reso- 
lution and sensitivity of the ALMA telescope, even in spectral ranges 
of low atmospheric transmission, are providing spatially and spec- 
trally resolved images of the vapour of the main water isotopologue 
in a planet-forming disk. Analyses of the morphology of the water 
emission, the spectrum of the highest SNR line and the excitation 
conditions jointly indicate that the (sub) millimetre lines are probing 
warm gas in the disk upper layers above the water snow surface, with 
aradially decreasing temperature profile. The non-detection of H80 
and the low water-to-dust ratio in the inner 17 au show that the obser- 
vations are probing marginally optically thick gas above the opaque 
dust continuum emission. These results highlight that the water 


content of quiescent protoplanetary disks at (sub) millimetre wave- 
lengths is most efficiently unveiled by targeting the main isotopo- 
logue, in particular in disks with high-continuum optical depths within 
the water snowline. 


Methods 

Observations, data reduction and imaging 

HL Tau was observed in both band 5 and band 7 with the ALMA 
program 2017.1.01178.S (principal investigator: E.H.), targeting the 
two para-water lines at 183.31004 and 325.15297 GHz, respectively. 
HL Tau was also observed in band 7 to target the latter water line 
with program 2022.1.00905.S (principal investigator: S.F.), together 
with the HO line at 322.46517 GHz (see Table 1). The molecular 
coefficients for the three transitions are reported in Tables 1 and 3. 

In the 2017.1.01178.S ALMA program, the source was observed in 
band 5 on 21 September 2018 for a total integration time of 46 min, 
with 43 antennas and baselines ranging between 15 mand 1.4 km. The 
weather conditions during the observation were exceptional, with a 
median precipitable water vapour (PWV) column during the observa- 
tions of -0.19 mm. JO423-0120 was used as amplitude and bandpass 
calibrator and JO510+1800 for phase referencing. The band 5 spectral 
setup had six spectral windows (SPWs), five of which targeted dif- 
ferent molecular transitions, including the 183 GHz water line, anda 
1.875-GHz-wide spw for continuum observations at 170.004 GHz. Inthe 
same program, HL Tau was observed in band 7 intwo spectral setups. 
The first was on 12 August 2019 with a time on-source of 31 min, with 
48 antennas and baselines ranging between 41 m and 3.6 km. The 
median PWV column was 0.4 mm.J0431+1731 was used to cross-calibrate 
phases, and JO538-4405 for amplitude and spectral response. This first 
spectral setup consisted of four SPWs in frequency division mode: three 
with a maximum bandwidth of 1.875 GHz, and one with 1,920 244 kHz 
channels, targeting the water 325 GHz line. The second spectral setup 
was used in the same program with an execution block observed on 
24 November 2017. The median PWV was 0.5 mm. JO519-4546 was 
used as amplitude and bandpass calibrator, and J0440+1437 as phase 
calibrator. The observation spent 33 min on the science target with 
49 antennas, with a maximum baseline of 8,500 m. The spectral setup 
consisted of four SPWs in frequency division mode: three with maxi- 
mum bandwidth, and one with 1,920 244 kHz channels, targeting the 
water 321 GHzline. 

Finally, new data were taken in October 2022 with a more compact 
array in band 7 with program 2022.1.00905.S, with baselines ranging 
between 15 and 500 m and a total on-source time of 100 min with 
41/45 antennas (in two execution blocks). The median PWV column was 
~0.3. J0423-0120 was used for flux and bandpass calibration, while 
J0431+1731 was used for phase referencing. The setup had four SPWs, 
one centred on the 325.15 GHz water line and another targeting the 
HO line at 322.46517 GHz. 

The cycle 6 (9) data were calibrated by the ALMA pipeline using 
CASA v.5.4 (v.6.4)". For the band 7 data, we combined the data from 
the two cycles. For both bands, we first self-calibrated each of the 
three execution blocks in phase, combining all spectral windows after 
flagging spectral ranges associated to line emission, and combiningall 
scans. Following ref. 43, we then aligned the data by fitting a Gaussian 
to the continuum and shifting the phase centre to the continuum 
maximum with the fixvis and fixplanet tasks. We then self-calibrated 
the two short-baseline (in the case of band 5, only one set of baselines 
is available) execution blocks in both phase and amplitude, reaching 
a peak SNR of 9,950 (a -500% improvement). In the case of band 7, we 
then combined the long-baseline execution block from the cycle 6 
data and self-calibrated the data again in both phase and amplitude, 
reaching a peak SNR of 5,300 (a 210% improvement). We took particular 
care with the amplitude calibration, where the gain solution 
with scan-length intervals greatly improved the data quality. 
The models for self-calibration were constructed with CLEAN with 
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Briggs weighting (robust = 0.5). The gain solutions were then applied 
to the full spectral data. 

The band 5 continuum data were imaged with robust = 0.0. With 
asynthesized beam of 0.364” x 0.312” (PA —9.5°), the band 5 (1.70 mm) 
continuum presents a r.m.s. of 33 wy, with a peak SNR of 2,457. The 
band 7 (0.94 mm) data were imaged with robust = -1.0. The data 
exhibit ar.m.s. of 36 y, with a peak SNR of 403 over asynthesized beam 
of 0.037” x 0.029” (PA 1.5°). A flux density for both images was obtained 
over an elliptical area with a semimajor axis of 1.3” and semiminor axis 
and PA to trace the disk inclination and PA (46.7 and 138.0°, respec- 
tively”). Without accounting for absolute flux calibration uncertain- 
ties, we obtain a flux density of 323.0 + 1.4 mJy and 1,677.7 + 0.6 mJy 
at 1.70 and 0.94 mm, respectively, where the uncertainties have 
been computed as standard deviations of randomly selected masks 
withthe same area of flux density extraction over emission-free regions 
of the continuum maps. 

While the water lines have been imaged with natural weighting (see 
main text), several additional attempts with a range of UV tapers were 
performed to increase the sensitivity to extended emission, but they did 
not show any feature undetected with natural weighting. Both integrated 
intensity (moment 0) and intensity-weighted velocity (moment 1) maps 
were generated for the water lines. Moment O maps were computed 
by integrating channels between -2 and 16 km s” without any clipping 
(Fig. 1) for the band 5 line, between -6.4 and 20.6 km s” for the 325 GHz 
line and between -10.4 and 24.6 km s” for the 321 GHz line. For the two 
brighter lines, we integrated the moment O maps over a circular area 
centred over the emission peak and a radius of 0.7”. We obtained line 
fluxes of 973 + 89 mJy kms“ and 929 + 89 mJy km s” (for the 183 and 
325 GHz lines, respectively). Since the 325 GHz line shows an absorp- 
tion red-shifted component, we also computed the underlying flux by 
considering the blue-shifted side only and multiplying it by a factor of 2. 
The resulting flux is 1,332 + 89 mJy kms“. The uncertainties on the line 
fluxes were computed by bootstrapping over 100 circular apertures in 
emission-free regions of the map and do not account for absolute flux 
calibration uncertainties. The same operation was applied tothe 321 GHz 
line using a smaller 0.28” radius extraction area. The resulting flux is 
679 +135 mJy kms“. The flux is much lower than the tentative detection 
by ref. 20 with the Submillimeter Array, where, however, the line is shifted 
by 30 km s” from the systemic velocity, indicating that the proposed 
emission may be originating from a large-scale flow that we filter out in 
our high-resolution data. For the 183 and 325 GHzlines, intensity-weighted 
velocity maps were generated using the bettermoments package** 
after applying 40 and 30 clipping to individual channels, respectively. 

The non-detection of the H}80 322 GHz line is shown in Extended 
Data Fig. 1. 


Fitting of the 183 GHz line spectrum 

Given the high SNR of the 183 GHz line, we fitted its spectrum by ana- 
lytically computing model spectra of a geometrically thin Keplerian 
disk, similarly to ref. 45. To do so, we assumed that the peak brightness 
temperature of the line decreases with radius as a power law: 

-q 


TR) = To( W 


D 
10au 


For a given radius R and azimuth @, we assumed that the line inten- 
sity follows a Gaussian distribution velocity: 


HMy(V — V y 
xp ( H K,proj ) , 


BT) a 
2kpT 


d2 


KR, $, v) (2) 


where B,(7) is the Planck function at temperature T, kg is the Boltzmann 
constant, d is the distance of HL Tau (140 pc), ym, is the mass of the 
water molecule and the Keplerian velocity projected along the line of 
sight Ux po can be written as follows: 


GM..\!2 
UK,proj = (=) sinicos @, (3) 


where iis the source inclination (46.7°) and G is the gravitational 
constant. We considered an optically thick limit when assuming a 
thermal broadening with a kinetic temperature equalling the bright- 
ness temperature. The flux density of the line can then be computed 
at every velocity v by integrating across the whole disk: 


21 pRou 
F(v) = cosi f [ KR, @, v)RARd@, (4) 
O Rin 


where cos iaccounts for the geometrical projection onthe sky. We fixed 
R to 0.1 au (but the model is not sensitive to this value for reasonably 
small radii) and sampled the disk with 150 points in radius and 550 in 
azimuth. We then convolved the models with a Gaussian kernel with the 
same channel width as in the data and sampled them at the same veloci- 
ties. We kept three free parameters in the fitting procedure: T}, q and 
Ru We fitted the spectrum shown in Fig. 2 with the emcee package 
using flat priors on the three free parameters: [10,1500] K, [0,3], 
[2,200] AU, respectively. We used 30 walkers, 1,000 steps of burn-in and 
1,000 additional steps to sample the posterior distribution. Extended 
Data Fig. 2 shows the best-fit model and 100 random draws extracted 
from the posterior distribution. While the fit does not manage to con- 
strain the outer radius of the emission, we obtain To = 287+!8° and 


—154 


q = 0.921020 The fit highlights a negative brightness temperature gradi- 
entin the radial profile, as seen in the reconstructed integrated intensity 


profile (Fig. 3) and as hinted by the rotational diagram shown in Fig. 4. 


Parametric fit of the integrated intensity profile in the 
visibility plane 

In order to extract the radial extent of the 183 GHz line, we performed 
a parametric fit of its integrated intensity radial profile in the visi- 
bility plane by exploiting the galario package. After averaging the 
continuum-subtracted visibilities in the same spectral range used to 
compute the moment O map, we fitted the visibility data by Fourier 
transforming a projected integrated intensity profile in the same uv 
points sampled during the observations. We modelled the radial profile 
with a simple Gaussian prescription 


R2 
JR) =Jo exp (-5). (5) 


G 


where we considered four free parameters: Jo, o, and the disk centre 
(ARA and ADec). We fixed the inclination and PA to those obtained from 
high-resolution continuum imaging”. The fit was performed with the 
emcee package, where the Jẹ parameter was sampled in log space. We 
used the following flat priors on the parameters: 


log,9U/steradian) € [1,20], 0, € [0,1.5]", ARA € 
[-0.4, 0.4], Adec. € [—0.4, 0.4]”. 


The posterior distribution was sampled with 50 walkers and 
1,000 steps, after 1,000 steps of burn-in. The MCMC exploration of 
the posterior space converges well, as shown in Extended Data Fig. 3. 


Rotational diagram and optical depth constraints 
To compute the rotational diagram of the water molecule, we use the 
same approach as in, for example, refs. 46,47. In the optically thin 
assumption, we can compute the column density Nnn and the rotational 
temperature 7,,, by measuring the integrated flux S,Av: 


N thin 


(6) 
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where Q is the solid angle used for flux extraction (see previous 
section), A,, is the Einstein coefficient of the considered transition, 
hand c are the Planck constant and the speed of light in vacuum, Q 
is the partition function, E, is the upper state energy (in K) and g, is 
the upper state degeneracy. Using the relation between N, tnin and Minin, 
the same equation can be written in logarithmic form’: 

In Nusnin 


=In Nihin — In QlTrot) — Eu/ Trot- (7) 


u 


We performed a linear regression using the emcee sampler” in 
the [InN(Nu tnin/Zu), Eu Space to extract N,nin and T,.,. We used the molecular 
coefficients reported in Table 3. In the fitter, the partition function was 
determined with a cubic spline interpolation across the rotational 
temperatures listed in ref. 37. The same approach was used for the 
rotation diagram analysis with six to eight transitions of o- and p-water 
in evolved stars*’. In the MCMC sampling, we used 128 walkers and 
2,000 steps (after 1,000 steps of burn-in). While Fig. 4 shows the result 
and the marginalized posterior distributions of the fit of all three lines, 
Extended Data Fig. 4 portrays the individual fits on the two colder and 
warmer lines, respectively. 

To compute the upper limit on the optical depth of the 325 GHz 
water line, we exploited the non-detection of the H150 line, which 
has almost identical molecular coefficients. Assuming that the column 
density of the main isotopologue line (t,o) is equal to 530 x TyI80 
(ref. 40), we can write 


(SyAV)1,0 1—e M0 
(SyAV)jpsq 1 — @=%,0/530" 


(8) 


Using the 3o upper limit on the HO transition and the measured 
flux of the H,O line (after correcting for absorption, since the absorbing 
column of H80 will have the same scaling factor of 530), a numerical 
solution of the equations leads to t,o < 14, as reported in the main 


text. 


Data availability 

All the ALMA data are publicly available on the ALMA archive (https:// 
almascience.nrao.edu/aq/), with program IDs 2017.1.01178.S and 
2022.1.00905.S. 


Code availability 

The Python packages used in the data analysis are all publicly 
available. The calibration and fitting scripts can be obtained from 
S.F. upon reasonable request. 


References 

1. Kreidberg, L. et al. A precise water abundance measurement for 
the hot Jupiter WASP-43b. Astrophys. J. Lett. 793, L27 (2014). 

2. Madhusudhan, N. Exoplanetary atmospheres: key insights, 
challenges, and prospects. Annu. Rev. Astron. Astrophys. 57, 
617-663 (2019). 

3. Rustamkulov, Z. et al. Early release science of the exoplanet 
WASP-39b with JWST NIRSpec PRISM. Nature 614, 659-663 
(2023). 

4. van Dishoeck, E. F., Bergin, E. A., Lis, D. C. & Lunine, J. I. in 
Protostars and Planets VI (eds Beuther, H. et al.) 835-858 
(Univ. Arizona Press, 2014). 

5. Drazkowska, J. et al. in Protostars and Planets VII (eds Inutsuka, S. 
et al.) 717-758 (Univ. Arizona Press, 2023). 

6. Cuzzi, J. N. & Zahnle, K. J. Material enhancement in protoplanetary 
nebulae by particle drift through evaporation fronts. Astrophys. J. 
614, 490-496 (2004). 

7. Schoonenberg, D. & Ormel, C. W. Planetesimal formation near the 
snowline: in or out? Astron. Astrophys. 602, A21 (2017). 


10. 


11. 


12. 


13. 


14. 


15. 


16. 


17. 


18. 


19. 


20. 


21. 


22. 


23. 


24. 


25. 


26. 


27. 


28. 


29. 


30. 


31. 


Drazkowska, J. & Alibert, Y. Planetesimal formation starts at the 
snow line. Astron. Astrophys. 608, A92 (2017). 

Oberg, K. l., Murray-Clay, R. & Bergin, E. A. The effects of 
snowlines on C/O in planetary atmospheres. Astrophys. J. 743, L16 
(2011). 

Eistrup, C., Walsh, C. & van Dishoeck, E. F. Setting the volatile 
composition of (exo)planet-building material. Does chemical 
evolution in disk midplanes matter? Astron. Astrophys. 595, A83 
(2016). 

Oberg, K. l., Facchini, S. & Anderson, D. E. Protoplanetary disk 
chemistry. Annu. Rev. Astron. Astrophys. 61, 287-328 (2023). 
Hogerheijde, M. R. et al. Detection of the water reservoir in a 
forming planetary system. Science 334, 338 (2011). 

van Dishoeck, E. F. et al. Water in star-forming regions: physics 
and chemistry from clouds to disks as probed by Herschel 
spectroscopy. Astron. Astrophys. 648, A24 (2021). 

Pontoppidan, K. M., Salyk, C., Blake, G. A. & Käufl, H. U. Spectrally 
resolved pure rotational lines of water in protoplanetary disks. 
Astrophys. J. Lett. 722, L173-L177 (2010). 

Grant, S. L. et al. MINDS. The detection of "CO, with JWST-MIRI 
indicates abundant CO, in a protoplanetary disk. Astrophys. J. Lett. 
947, L6 (2023). 

Kóspál, Á. et al. JWST/MIRI spectroscopy of the disk of the young 
eruptive star EX Lup in quiescence. Astrophys. J. Lett. 945, L7 
(2023). 

Banzatti, A. et al. JWST reveals excess cool water near the snow 
line in compact disks, consistent with pebble drift. Astrophys. J. Lett. 
957, L22 (2023). 

Salyk, C. et al. Detection of water vapor in the terrestrial planet 
forming region of a transition disk. Astrophys. J. Lett. 810, L24 (2015). 
Tobin, J. J. et al. Deuterium-enriched water ties planet-forming 
disks to comets and protostars. Nature 615, 227-230 (2023). 
Kristensen, L. E., Brown, J. M., Wilner, D. & Salyk, C. Velocity- 
resolved hot water emission detected toward HL Tau with the 
Submillimeter Array. Astrophys. J. Lett. 822, L20 (2016). 

Carr, J. S., Najita, J. R. & Salyk, C. Measuring the water snow line in 
a protoplanetary disk. Res. Notes AAS 2, 169 (2018). 

Bosman, A. D. & Bergin, E. A. Reimagining the water snowline. 
Astrophys. J. Lett. 918, L10 (2021). 

Notsu, S. et al. Dust continuum emission and the upper limit fluxes 
of submillimeter water lines of the protoplanetary disk around HD 
163296 observed by ALMA. Astrophys. J. 875, 96 (2019). 

ALMA Partnership. et al. The 2014 ALMA long baseline campaign: 
first results from high angular resolution observations toward the 
HL Tau region. Astrophys. J. Lett. 808, L3 (2015). 

Yen, H.-W. et al. HL Tau disk in HCO* (3-2) and (1-0) with ALMA: gas 
density, temperature, gap, and one-arm spiral. Astrophys. J. 880, 
69 (2019). 

Riviere-Marichalar, P. et al. Detection of warm water vapour in 
Taurus protoplanetary discs by Herschel. Astron. Astrophys. 538, 
L3 (2012). 

Salyk, C. et al. A high-resolution mid-infrared survey of water 
emission from protoplanetary disks. Astrophys. J. 874, 24 (2019). 
Belitsky, V. et al. ALMA band 5 receiver cartridge. Design, 
performance, and commissioning. Astron. Astrophys. 611, A98 
(2018). 

Notsu, S. et al. Candidate water vapor lines to locate the H,O 
snowline through high-dispersion spectroscopic observations. Il. 
The case of a Herbig Ae star. Astrophys. J. 836, 118 (2017). 

Notsu, S. et al. Candidate water vapor lines to locate the HO 
snowline through high-dispersion spectroscopic observations. III. 
Submillimeter H,O and H,O lines. Astrophys. J. 855, 62 (2018). 
CASA Team. et al. CASA, the Common Astronomy Software 
Applications for radio astronomy. Publ. Astron. Soc. Pac. 134, 
114501 (2022). 


Nature Astronomy 


Article 


https://doi.org/10.1038/s41550-024-02207-w 


32. Garufi, A. et al. ALMA chemical survey of disk-outflow sources in 
Taurus (ALMA-DOT). V. Sample, overview, and demography of disk 
molecular emission. Astron. Astrophys. 645, A145 (2021). 

33. Garufi, A. et al. ALMA chemical survey of disk-outflow sources in 
Taurus (ALMA-DOT). VI. Accretion shocks in the disk of DG Tau 
and HL Tau. Astron. Astrophys. 658, A104 (2022). 

34. Tazzari, M., Beaujean, F. & Testi, L. GALARIO: a GPU accelerated 
library for analysing radio interferometer observations. Mon. Not. 
R. Astron. Soc. 476, 4527-4542 (2018). 

35. Teague, R. GoFish: fishing for line observations in protoplanetary 
disks. J. Open Source Soft. 4, 1632 (2019). 

36. Carrasco-Gonzalez, C. et al. The radial distribution of dust 
particles in the HL Tau disk from ALMA and VLA observations. 
Astrophys. J. 883, 71 (2019). 

37. Polyansky, O. L. et al. ExoMol molecular line lists XXX: a complete 
high-accuracy line list for water. Mon. Not. R. Astron. Soc. 480, 
2597-2608 (2018). 

38. Bethell, T. & Bergin, E. Formation and survival of water vapor in the 
terrestrial planet-forming region. Science 326, 1675 (2009). 

39. Baulch, D. L. Evaluated Kinetic Data for High Temperature 
Reactions (CRC Press, 1972). 

40. McKeegan, K. D. et al. The oxygen isotopic composition of the sun 
inferred from captured solar wind. Science 332, 1528 (2011). 

41. Schroeder I, I. R. H. G. et al. °O/"°O ratio in water in the coma of 
comet 67P/Churyumov-Gerasimenko measured with the Rosetta/ 
ROSINA double-focusing mass spectrometer. Astron. Astrophys. 
630, A29 (2019). 

42. Weaver, E., Isella, A. & Boehler, Y. Empirical temperature measure- 
ment in protoplanetary disks. Astrophys. J. 853, 113 (2018). 

43. Andrews, S. M. et al. The Disk Substructures at High Angular 
Resolution Project (DSHARP). I. Motivation, sample, calibration, 
and overview. Astrophys. J. 869, L41 (2018). 

44. Teague, R. & Foreman-Mackey, D. A robust method to measure 
centroids of spectral lines. Res. Notes AAS 2, 173 (2018). 

45. Zagaria, F. et al. Testing protoplanetary disc evolution with 
CO fluxes. A proof of concept in Lupus and Upper Sco. 

Astron. Astrophys. 672, L15 (2023). 

46. Loomis, R. A. et al. The distribution and excitation of CH,CN ina 
solar nebula analog. Astrophys. J. 859, 131 (2018). 

47. Facchini, S. et al. The chemical inventory of the planet-hosting 
disk PDS 70. Astron. J. 162, 99 (2021). 

48. Goldsmith, P. F. & Langer, W. D. Population diagram analysis of 
molecular line emission. Astrophys. J. 517, 209-225 (1999). 

49. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: 
the MCMC hammer. Publ. Astron. Soc. Pac. 125, 306 (2013). 

50. Baudry, A. et al. ATOMIUM: probing the inner wind of evolved 
O-rich stars with new, highly excited H,O and OH lines. 

Astron. Astrophys. 674, A125 (2023). 

51. Pickett, H. M. et al. Submillimeter, millimeter and microwave 
spectral line catalog. J. Quant. Spectrosc. Radiat. Transf. 60, 
883-890 (1998). 

52. Barber, R. J., Tennyson, J., Harris, G. J. & Tolchenov, R. N. A 
high-accuracy computed water line list. Mon. Not. R. Astron. Soc. 
368, 1087-1094 (2006). 


Acknowledgements 

This paper makes use of the following ALMA data: ADS/JAO. 
ALMA#20171.01178.S, ADS/JAO.ALMA#2022.1.00905.S. ALMA is a 
partnership of the European Southern Observatory (representing 
its member states), the National Science Foundation (USA) and 
the National Institutes of Natural Sciences (Japan), together with 
the National Research Council (Canada), the Ministry of Science 
and Technology and Academia Sinica Institute of Astronomy and 


Astrophysics (Taiwan), and the Korea Astronomy and Space Science 
Institute (Republic of Korea), in cooperation with the Republic of 
Chile. The Joint ALMA Observatory is operated by the European 
Southern Observatory, Associated Universities Inc./National Radio 
Astronomy Observatory and the National Astronomical Observatory 
of Japan. S.F. is funded by the European Union (European Research 
Council (ERC), UNVEIL, 101076613). Views and opinions expressed, 
however, are those of the author(s) only and do not necessarily 
reflect those of the European Union or the ERC. Neither the European 
Union nor the granting authority can be held responsible for them. 
S.F. acknowledges financial contribution from PRIN-MUR POPS 
2022YP5ACE. L.T. acknowledges funding from Progetti Premiali 2012 
iALMA (CUP C52I13000140001), Deutsche Forschungs-gemeinschaft 
(German Research Foundation) ref no. 325594231 FOR 2634/1 TE 
1024/1-1, European Union’s Horizon 2020 research and innovation 
programme under the Marie Sklodowska-Curie grant no. 823823 
(DUSTBUSTERS) and the ERC via the ERC Synergy Grant ECOGAL 
(grant no. 855130). M.V.D. acknowledges support from Wallonie- 
Bruxelles International (Belgium) through its grant ‘Stage en 
Organisation Internationale’ and the French-speaking Community 
of Belgium through its FRIA grant. 


Author contributions 

S.F. and E.H. led the ALMA observing proposals. S.F. calibrated the 
data and performed the data analysis. S.F. and LT. wrote the manuscript. 
All co-authors provided input on the manuscript. 


Competing interests 
The authors declare no competing interests. 


Additional information 
Extended data is available for this paper at 
https://doi.org/10.1038/s41550-024-02207-w. 


Correspondence and requests for materials should be addressed to 
Stefano Facchini. 


Peer review information Nature Astronomy thanks Shota Notsu and 
the other, anonymous, reviewer(s) for their contribution to the peer 
review of this work. 


Reprints and permissions information is available at 
www.nature.com/reprints. 


Publisher's note Springer Nature remains neutral with regard 
to jurisdictional claims in published maps and institutional 
affiliations. 


Open Access This article is licensed under a Creative Commons 
Attribution 4.0 International License, which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, 
as long as you give appropriate credit to the original author(s) and the 
source, provide a link to the Creative Commons licence, and indicate 
if changes were made. The images or other third party material in this 
article are included in the article’s Creative Commons licence, unless 
indicated otherwise in a credit line to the material. If material is not 
included in the article's Creative Commons licence and your intended 
use is not permitted by statutory regulation or exceeds the permitted 
use, you will need to obtain permission directly from the copyright 
holder. To view a copy of this licence, visit http://creativecommons. 
org/licenses/by/4.0/. 


© The Author(s) 2024 


Nature Astronomy 


Article https://doi.org/10.1038/s41550-024-02207-w 


— 1 km s™! channel 
—— 10 km s7! channel 


=== Vsys 


Flux Density (J y) 


-20 -10 0 10 20 30 
Velocity (km/s) 


Extended Data Fig. 1 | Non-detection of the 322 GHz H380 line. Spectrum of the 322 GHz H80 line, extracted over a circle with radius of 0. 7’ centered on the 


continuum peak, as for Fig. 2. Two spectra are shown, from channel maps with two different channel widths. The grey dashed line indicates the systemic velocity of 
7.1kms*, In the top left 2oscale bars are reported for reference. 
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Extended Data Fig. 2 | Corner plot of the radial extent of the 183 GHz line. Marginalized posterior distribution of the galario fit of the visibility points of the 
integrated intensity of the 183 GHz line (see Fig. 3). The dashed lines indicate the 16th, 50th and 84th percentiles of the marginalized posterior distributions. 
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Extended Data Fig. 3 | Fit of the 183 GHz line spectrum. Left: spectrum of the the dark red line indicating the bestfit model. Right: marginalized posterior 
183 GHz H,O line, as in Fig. 2. The uncertainty associated to each data point is distribution of the fitted parameters, where no constraint can be retrieved 
shown in the left part of the spectrum as an errorbar. A random sampling of on Rou from this analysis. 
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Rotational diagram for two colder lines 
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Rotational diagram for two warmer lines 
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Extended Data Fig. 4 | Rotational diagrams of line pairs. Rotational diagrams dashed line in the left panels shows the best fit. The black dashed lines in the 
and marginalized posterior distribution of the MCMC exploration for individual corner plots indicate the 16th, 50th and 84th percentiles of the marginalized 
fits on the two lowest and two highest energy lines, respectively. The green posterior distributions. 


Nature Astronomy 


